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Abstract 

This paper presents a strategy for a posteriori error estimation for substruc- 
tured problems solved by non-overlapping domain decomposition methods. 
We focus on global estimates of the discretization error obtained through the 
error in constitutive relation for linear mechanical problems. Our method al- 
lows to compute error estimate in a fully parallel way for both primal (BDD) 
and dual (FETI) approaches of non-overlapping domain decomposition what- 
ever the state (converged or not) of the associated iterative solver. Results 
obtained on an academic problem show that the strategy we propose is effi- 
cient in the sense that correct estimation is obtained with fully parallel com- 
putations; they also indicate that the estimation of the discretization error 
reaches sufficient precision in very few iterations of the domain decomposi- 
tion solver, which enables to consider highly effective adaptive computational 
strategies. 

Key words: verification; error in constitutive relation; non overlapping 
domain decomposition; FETI; BDD. 



1. Introduction 

The setting-up of robust numerical methods to solve complex systems of 
partial differential equations has become a key issue in applied mathematics 
and engineering, driven by the increasing use of numerical simulation in both 
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research and industry. Among the latter, virtual testing has become a short 
term aim, with the objective to replace expensive experimental studies and 
validations by numerical simulations, even in order to certify large structures 
as planes and bridges. 

Thus, one key point of the numerical methods to develop is the verification 
of computations which enables to warranty that the computed solution is 
sufficiently close to the original continuum mechanics model. This topic of 
numerical analysis has been the subject of many studies for the last decades. 
Three main classes of error estimator have been developed, based either on 
equilibrium residuals [H, flux projection j2| or error in constitutive law [3[. 
An overview of those various methods can be found in |4(. 

Another key point of numerical methods is their ability to quickly provide 
solutions to large (nonlinear) systems. The most classical answer to this issue 
is to use domain decomposition methods in order to take advantage of the 
parallel hardware architecture of recent clusters and grids. In engineering, 
non-overlapping domain decomposition methods are mostly employed, such 
as the well known FETI or BDD [6] . An overview of the main approaches 
related to non-overlapping domain decomposition can be found in [7J. 

We aim to provide fully integrated adaptive strategies to compute large 
structural mechanics problems with certified quality. To do that, our cur- 
rent approach is to explore some ways of making bidirectional interactions 
between domain decomposition and a posteriori error estimation. Our de- 
velopments are based both on the error in constitutive relation to measure 
the quality of our results and to forecast mesh refinement, and on a generic 
vision of non-overlapping domain decomposition methods which enables to 
do high-performance computing. 

This paper focuses on the estimation of the global error in constitutive 
relation in order (among others) to study how it is influenced by the error 
in the convergence of the domain decomposition solver which is linked to 
the non-satisfaction of interface equations (continuity of displacements and 
balance of forces). To do so we propose a strategy to build, in parallel and 
during the iterations, displacement and stress fields which are kinematically 
admissible (KA) and statically admissible (SA) on the whole structure. We 
face two main difficulties. First, since before convergence interface fields do 
not possess the classical properties of discretized fields (continuity of displace- 
ments and weak equilibrium), the recovery of admissible displacements and 
stresses requires some preprocessing. Second, the computation of statically 
admissible fields being an operation which can not be conducted indepen- 
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dently on each element (in some methods it can even be a large bandwidth 
operation), classical recovery methods jil, 0, [lo[ U, 12, 13] would require 
inter-subdomain communications. 

Our generic method to build continuous displacement and balanced trac- 
tion fields for both primal and dual approaches of non-overlapping domain 
decomposition is presented through this paper. It will be shown that the 
properties of the preconditioners involved in domain decomposition solvers 
make this reconstruction costless, and that an error estimator can then be 
computed in a fully parallel way. 

This paper is organized as follows. Section [2]recalls the general framework 
related to our upcoming developments, mainly the estimation of the error in 
constitutive equation and the use of domain decomposition method. Section 
|3] shows how the problem of error estimation in a substructured context can 
be brought back to the computation of nodal displacement and traction fields 
which are admissible in a discrete sense. Sections H] and |5] describes how to 
obtain these fields without inter-sub domains exchanges when using classical 
primal (BDD) and dual (FETI) domain decomposition methods with good 
preconditioners. Section [H] presents numerical assessments, first to validate 
the parallel recovery procedure, then to prove that a good estimation can 
be obtained far earlier than the solver converged (in the sense of domain 
decomposition iterative solver). Finally, Section [7] concludes this paper. 



2. Framework of the study 

2.1. Reference mechanical problem 

Let us consider the static equilibrium of a structure which occupies the 
open domain Q C M. d and which is submitted to given body forces /, to given 
traction forces g on dffl and to given displacements uq on the complementary 
part d u Q ^ 0. We assume the structure undergoes small perturbations and 
that the material is linear elastic, characterized by the Hooke's tensor EL 
Let u be the unknown displacement field, e(u) the symmetric part of the 
gradient, a the Cauchy stress tensor. 

Let lo C fl be an open subset of fl, dfUJ = dui fl dffl, d u u = du fl d u fl and 
r = du \ (d u u U dfu) (see Figured]). We introduce two affine subspaces and 
one positive form: 

• Subspace of kinematically admissible fields 

KA(w) = {ue (HV)) d , tr(u) = u on d u u) (1) 
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Figure 1: Domain fi, subdomain u) and boundaries 

where tr is the trace operator. 
Subspace of statically admissible fields 

SA(w) = |r G (L 2 (u)) dxd ,r symmetric, Wu* G KA°V) , 

/ r : e{u*)duj = / f.u*du + / g.u*dS 

Juj Jul JdfU) 

where KA 00 (u;) = |w G (H 1 ^)) 01 , tr(«) = on d u w U r} 
Measure of the non-verification of the constitutive equation 
ecR(u/)0,o") = ||(T-1HI: e(u) \\w-\u 



where llxlL-i u , = \ / (x : HI -1 : a;) du 



The mechanical problem set on Q can be formulated as: 

Find (u ex ,a ex ) G KA(fi) x SA(fi) such that e C R(n)(w ea; , o ex ) = 



2.2. Finite element approximation for the global problem 

Let Qh be a tessellation of Cl to which we associate a finite dimensional 
subspace KA^(f2) of KA(fi). The classical finite element displacement ap- 
proximation consists in searching 

u h G KA h (n) 
a h = H : e{u h ) 

a h : e(u* h )dQ = [ f.u* h dQ + [ g.u* h dS, V< G KA° h °(Q) 



n Jn Jd f n 

After introducing the d x A^o/ matrix ip h of shape functions which form 
a basis of KA/j(f2) and the vector of nodal unknowns u (of size A^o/, number 
of degrees of freedom) so that Uh = ¥>h u > the classical finite element method 
leads to the well-known linear system: 

Ku = f (5) 

where K is the (symmetric positive definite) stiffness matrix of domain fl^ 
and f is the vector of generalized forces. 

2.3. A posteriori error estimator 

The finite element approximation (uh, ah) satisfies Uh G KA(Q) and 
e CR(a)(uh, o~h) = but ah ^ SA(fi). The error in constitutive relation 
consists in deducing from (uh, ah) an admissible displacement-stress pair 
(uh,ah) G KA(Q) x SA(fi) in order to measure the residual on the consti- 
tutive equation ([3]) ecR(o)(«/i, a h ) ^ 0. Using the well-known Prager-Synge 
theorem it can be proved that 

\\e(u ex ) -e(wfc)llH,n+ Ike* - ^lln-i.n = { e CR(n){u h ,Vh)) 2 

Hence, the evaluation of the error in constitutive relation ecR(n)(uh,&h) for 
any admissible pair (uh, a~h) provides a guaranteed upper bound of the global 
error 

\\e{u ex ) -e(u A )||H,n < e CR{n) (uh,a h ) (6) 

KA/ l (f2) being a subspace of KA(f2), the construction of an admissible dis- 
placement field Uh is straightforward since it can be taken equal to Uh- On the 
other hand, as ah is not statically admissible, the construction of an admis- 
sible stress field dh G SA(f2) is a crucial point which has already been widely 
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studied in the literature. A first solution is to use a dual formulation of the 
reference problem [l4| to compute a from scratch. Unfortunately building a 
subspace of SA(f2) is a complex task and most people prefer to post-process 
a statically admissible field from Field cr^ obtained by a displacement for- 
mulation. Classical methods are the element equilibration techniques 
which have been improved by the use of the concept of partition of unity 
which lead to [HI, [HI, [l3| and the flux-free method jioj]. In most cases they 
involve the computation of efforts on "star-patches" which are the set of el- 
ements sharing one node, for each node of the mesh. Though rather simple 
these computations are in great number and thus expensive. 

In the following, we note by J-'h the algorithm which has been chosen to 
build an admissible stress field <5v Whatever the choice, the algorithm takes 
as input not only the finite element stress field ah but also the continuous 
representation of the imposed forces (f,g)- 

bh = Fhivh, f, g) e SA(fi) 

The algorithm we have used for our applications is the one proposed in ^ 
using a three deg rees higher polynomial basis when solving the local problems 



on elements 15]. 



2.4- Substructured formulation 

Let us consider a decomposition of domain fl in open subsets {^ s ')i^ s ^N Bd 
(N s d is the number of subdomains) so that ft^ s > fl Cl^ s -* = for s ^ s' and 
Q = U S Q( S \ Let u u = (u^) s , we define the global assembling operator A: 

u = A(u n ) <^> U| nW = u {s) (7) 

In order to reformulate the mechanical problem on the substructured con- 
figuration, we need to specify the conditions that should be satisfied at the 
boundary between subdomains r (ss '- ) = dfl^ fl dil^ s '\ We have the funda- 
mental properties: 



A( nx J uis) e KA(fiW) , Vs 

' E KA ^ j ^ \ tr(«W) = tr(u(*')) on T^'\ V(s, s') 

At nx J a(s) G SA(fi^) , Vs 

A{a ) E bA{U) O | ^ w + n(0 = Q Qn p K)) v(s? ^ 



(8) 
(9) 
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In other words, in order to be admissible on the whole domain fi, not only 
fields need to be admissible in a local sense (independently on each f2^), but 
they also need to satisfy interface conditions, namely displacements continu- 
ity and tractions balance (action-reaction principle). 

2.5. Finite element approximation for the sub structured problem 

We assume that the tessellation of Cl and the substructuring are conform- 
ing so that (i) each element only belongs to one subdomain and (ii) nodes are 
matching on the interfaces. Each degree of freedom is either located inside 
a subdomain (subscript i) or on its boundary = U S >T( SS > (subscript b) 
where it is shared with at least one neighboring subdomain. Let A^ be the 
vector of unknown efforts imposed on the interface of subdomain fl^ by its 
neighbors. The finite element problem (jSJ) can be written highlighting the 
contributions of subdomains: 



Ea^a^ = o 



Ws, K^uW = f« + t« J A« with £ a(s)u ( s) = q (10) 

V s 

where is the discrete trace operator (u£ = t^u^) and where 
and are assembling operators so that A^ s ^ enables to formulate the 
mechanical equilibrium of interfaces (Q and A^ enables to formulate the 
continuity of displacements (jSJ (in the case of two subdomains, we have 
J2 S A (s) X b = A[ 1} + A[ 2) = and £ s A«u< a) = - u[ 2) = 0, see Fig. [flfor 
less trivial example and |7| for more an extensive description of all operators). 
One fundamental property of assembling operators is their orthogonality: 

^a«am t = o (11) 

s 

Note that the equilibrium of subdomain also writes: 

(12) 

7 \H J \ A b y 

or in an equivalent condensed form: 

S^u[ s) = + a[ s) (13) 
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with 

where SW is the Schur complement and b( s ) is the condensed right-hand side. 



3. A posteriori error estimator in substructured context 

The key point for the efficient evaluation of the error in constitutive rela- 
tion in a substructured context (without overlapping) is to define admissible 
pairs (u^\a^) € KA(f^) x SA(f^) on each subdomain so that the associ- 
ated assembled pair is admissible for the reference problem (A(v^), A(a^)) G 
KA(Q) x SA(fi). Due to the absence of overlap, the additive structure of the 
associated error in constitutive relation leads to a fully parallel evaluation of 
the a posteriori error estimator: 

(e C R(n){A(^),A(dh))) 2 = Yl ( e CR(f^)) ^ s) )) 

s 

The application of a classical recovery strategy to compute admissible 
fields raises two difficulties in a substructured context. First, the star-patches 
can not be employed on the boundary nodes without assuming communica- 
tion between subdomains. Though these exchanges would remain limited, 
we propose an alternate strategy to achieve full parallelism without impair- 
ing the properties of the error in constitutive relation. Second, in order to 
solve the substructured problem (fit)]) parallel strategies consist in using iter- 
ative solvers which are based on the loosening of at least one of the interface 
conditions which is only verified (up to a certain precision) once the solver 
converged. Thus recovering strategies need to be adapted so that the local 
fields (uj*\dfr) satisfy the interface conditions. 

The aim of this section is to prove that the determination of the admissi- 
ble pair (A(u^), A(&h)) can be brought back to the determination of nodal 

interface fields (u£ s \A^) s which satisfy specific interface conditions. The 
construction of these nodal fields depends on the chosen domain decomposi- 
tion strategy and is discussed in the following sections. 
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3.1. Kinematically admissible fields 

In order to ensure interface Condition (JS} when building G KA(f2^) 
so that A(v^) G KA(f2), we introduce continuous interface displacement 

A (S) 

fields from which we shall deduce internal displacement fields: 



= fig? , v(«y) 

?y (s) - ?y (s) Vs 
fc|r( ss ') ~~ 6^ ' 

Since discretizations are matching on the interface, the first condition can 
directly be imposed on finite element nodal quantities: 

u< s) = uj s '\ V( S , S ') 

In order to deduce the internal fields, one finite element problem is solved 
independently on each subdomain with imposed Dirichlet conditions on the 
interface: 



a!" = kF 



U 6 



tt = -4(^) G KA(fi) 

S 1 .^. Statically admissible fields 

In order to ensure interface Condition ([5]) when building a^' G SA(f2^) 
so that A(Oh) G SA(O), we introduce for each subdomain the continuous 
balanced interface traction fields defined on which satisfy: 

^ ] -n {s) = on rW 

*£> + = o on r(-') (15) 

/ f.pdtt+ [ g('\pdS+ [ F$.pdS = Vp G RKA°(fiW) 

where RKA°(f2( s )) is the set of rigid body motions which are compatible with 
Dirichlet conditions imposed on d u Cl^: 

RKA°(fi (s) ) = {pe H^fiW), p = on d u ^ s) , e(p) = 0, } 
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The last condition of (I15p is the translation of Fredholm's alternative in 
order to ensure the well-posedness of the static problem on domain fiw. To 
build these traction fields in a simple way, we associate them with the finite 



element nodal reaction field X b : 

/ ,^-^\ r ^dS = \% (16) 
where j denote a node of the interface, ipj its associated shape function and 

~(s) ~(s) 

\ b j the corresponding nodal component of X b . This equation then imposes 

-~(s) ^(s) 

that the discrete field X b and the continuous field F bh develop the same 
virtual work in any finite element displacement field. The conditions on F bh 

have these discrete counterparts on X b : 

£ A( ^ S) = ° 

(17) 

R (s) T (t(°) T \ ( b s) + f M) = o 

where is a basis of ker(K^^). As said earlier, the first equation corre- 
sponds to the equilibrium between subdomains. The second equation corre- 
sponds to the balance of the subdomain with respect to virtual rigid body 
motions (since this kind of displacement field is exactly represented in the 
finite element approximation, the discrete condition is equivalent to the con- 
tinuous one). 

As a first approach, we define F b ^ as: 

F« = ,F« (18) 

"^(b) ~^(s) (s) 

where F b is the vector of nodal values of F bh and <Ph) T [ a ) refers to the vector 

of the trace on of finite element shape functions. Vector F b is then 
obtained by the inversion of the (small) "mass" matrix of the interface of 
each subdomain. In the following, we denote by Qh the previous procedure 
which associates a continuous balanced interface force F b ^ to a balanced 

~(s) 

nodal interfaces forces X b : 

F b (s h ] = e h & s) ) (is) 
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The traction field F bh ^ allows to satisfy the interface conditions associated 
to the static admissibility. The next step is to build internal finite element 

~(s) 

stress fields which match the associated nodal boundary field X b . This is 
done by solving one finite element problem on each subdomain with imposed 
Neumann conditions on the interface. 

u« = K« + (f W + t^A^) 

^=J-,(H: £ (^uW),/«,{^),^(Ah}) ^ 
a h = Affi) e SA(fl) 

The use of the pseudo-inverse is due to the potential lack of Dirichlet 

boundary conditions on the substructure. Displacement field is defined 
up to a rigid body motion which needs not to be determined since only the 
symmetric gradient of the associated displacement field is required. 

It has to be noted that the fully parallel procedure Gh proposed above 
leads to a different admissible traction field as would have been obtained 
using standard patch-technique ^ (referred in the sequel as the sequential 
approach). Thus the use of Gh implies that the parallel error estimation 
is different from the standard sequential one even when discrete interface 
conditions are satisfied. For now there are no theoretical results on the 
quality of the resulting fields, examples (as given in Section [6]) show that 
sequential estimator and parallel estimator (when interface conditions have 
sufficiently converged, which happens very quickly) can not be distinguished. 

4. Recovery of admissible fields in BDD 



In the Balancing domain decomposition [16|, Il7| . a unique interface dis- 



placement unknown U(, is introduced so that continuity is always insured: 
u[ s) = A^ T u b =► A (s) ui s) = (21) 

s 

Other quantities can be deduced from u b and equations (I12II13P : 
\[ s) = S^u b - b« 



(22) 
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The BDD solver consists in iteratively finding the interface displacement U;, 
which insure global equilibrium (Y^ s A^A^ = 0), 



o = A(s)x l s) = (x> (s)s(s)A(s)T ) u * - ( E A(s)b(s) 

S \ S / \ S / 



(23) 



4-1. Recovery of KA fields 

In the BDD solver, kinematic interface conditions are satisfied anytime 
and using = u;, enables to build uf 3 so that Uh = A (v^) e KA(fi). Note 
that all associated computations are realized during the standard resolution 
process so that no extra operation is required. 

4-2. Recovery of SA fields 

For a given interface displacement u b , we note: 

M = £ a«a<'> = A(s) ( s(s) ^ - b(s) ) 

s s 

Obviously |JAfoJ| is zero if and only if u b is the solution to (123"]). We then 
define: 

A^A^- A« T M (24) 

where (A^) s are scaled assembling operators so that ^ s A^) A^) = I. 
The multiplicity scaling is a typical example of such operator AW; 



~ t 1 

which, in the case of two subdomains, gives A^ [fA&]] = — ff Af,7l • I* 1 t ne 
case of heterogeneous structures, other scaled assembly operators which take 



the heterogeneity into account are used [18|, [19|, [20 

~(s) 

It is clear that by definition, X b is a balanced nodal reaction field: 

s 

~(s) 

In order to prove that \ b also satisfies Fredholm's alternative, we note that 
since is a basis of ker(K^) and is invertible, we have S^R^ = 
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and R^ = — K.\f K-^R^. The condition then writes in an equivalent 
condensed form: 

Rr T (Al" + b"))=0 

R<' ,T (s w u t - b<'> + A'" T IfAftH + b«) = 

Using the symmetry of (inherited from the symmetry of K^) to nullify 
S^ sj , the condition writes: 

(A^Ri s) ) T M = ( 25 ) 



which is exactly the balancing condition [6] of the iterative BDD solver: the 
residual of the BDD iterative solver f[A 6 "fl = fe s A^A[ s) ) ([23]) has to be 
orthogonal to all local weighted rigid body motions so that preconditioning 
step is well posed. 

Then we have constructed a pair of interface nodal Vectors (u&, A&) which 
satisfy all required conditions to build admissible fields. 

Note that all the involved operations are already realized during clas- 
sical steps of the primal domain decomposition approach with a Neumann- 
Neumann preconditioner and the associated coarse problem, so that all finite 
element quantities (even the internal ones) are available at no cost; the only 
extra operations are due to the use of Algorithms Qh (to compute F bh ) and 
Th (to compute ah). 

5. Recovery of admissible fields in FETI 

In the Finite Element Tearing and Interconnecting domain decomposi- 
tion [5] , a unique interface effort unknown A& is introduced so that interface 
equilibrium is always insured: 

\ { b s) = A« T A b A(s)A i S) = ( 26 ) 

s 

Displacements can be deduced from X b if it satisfies Fredholm's alternative 
on each substructure: 

U W = (f« + t^ T A^ T \ b ) + R^aW 

V (27) 
= RW T (f « + t^ T A^ T \ b ) 
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where ct^ is the unknown magnitude of rigid body motions. The FETI 
solver consists in iteratively finding an interface effort A&, under the previous 
constraint, which insures the continuity of interface displacement: 

o = J2 A (s)u i s) = fe a^k^\^ t a^A \ b 

+ ^A^K^A + ^AW'WwA (28) 

5.1. Recovery of SA fields 

In the FETI solver, the nodal interface fields A^ = A^ T A& are by 
construction always balanced at the interface (1261) and associated to well- 
posed discrete Neumann problems on each substructure (l2Tj) . Hence, we 



can directly set a[ ^ = A^ and apply algorithms Gh an d J~h to compute 
d ( h s) e SA(Q^) with a h = Affi) e SA(fl). 

5.2. Recovery of KA fields 

For a given balanced nodal interface traction A&, we introduce, in agree- 
ment with ( 1271) . the gap of the interface displacement : 

s 

and we define 

~( s ) 00 aW T m II 

u b = u l - A M 

~ (s) / \ ~ (s) T 

where (A ) a are scaled assembling operators so that Y]„ A^ s) A = I. 

~ (s) 

Similarly to the BDD typical example of such operator A is the 

multiplicity scaling: 

A {S)T = aW^TA^ A U)T ) (29) 



j 



~ (s) T 1 

Note that in the case of two subdomains, we have: A [[ubJJ = — [|u&J|- 
The connection between FETI and BDD scaling operators (even in the het- 



erogeneous case) is given in 21 
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It is clear that by construction 

s 

Hence nodal interface displacement can be used to deduce an admissible 
displacement field u h so that Uh = A (u^) G KA(fi). 

Then we have constructed a pair of interface nodal Vectors (uf,, A&) which 
satisfies all required conditions to build admissible fields. Note that all the 
involved operations are already realized during classical steps of the dual do- 
main decomposition approach (with built-in coarse problem) with Diri chief's 
preconditioner, so that all finite element quantities (even the internal ones) 
are available at no cost: the quantity [[u^JJ is directly available during the 
classical solution procedure (without computing any afity which is based on 
an initialization/projection algorithm (5)], and the displacement field u^ s ^ can 
be defined up to an element of the kernel (a rigid body motion) since only 
its symmetric gradient is used during the computation of the error. The only 
extra operations are due to the use of algorithms Gh (to compute Fbh) and 
Th (to compute ah). 



6. Numerical assessment 



In order to assess the performance of our parallel error estimator, we 



consider the 2D toy problem of the T-shape structure of Figure 2(a) which has 
been used in other papers like 22] . Plane stresses are assumed. The material 
behavior is isotropic, linear and elastic, with Young modulus E = 2000 MPa 
and Poisson's ratio v = 0.3. The structure is clamped on its basis (whose 
length is denoted L) and it is submitted to traction and shear on its upper- 
right side, while all the remaining boundaries are traction-free. 

Several regular meshes have been generated, constituted by triangular 
elements of characteristic size h = - with m = 2,4,8,16,32. For each 

m 

mesh, a sequential (mono-domain) computation is driven, followed by domain 
decomposition computations obtained by an automatic splitting of the mesh 
in an increasing number N s d of subdomains ( N s d = 2, 4, 8 when m ^ 4 and 
N sd = 2,4,8, 16,32 when m ^ 8). Figure 2(b) shows such a decomposition 
for N s d = 8 and m — 8. 



All the computations are driven in the ZeBuLoN finite element code [23 
using elements of polynomial degree p = 1. Both BDD and FETI algorithms 



15 



(a) Finite element problem (b) Substructuring 

(h= §) (h= ±,N sd = 8) 



Figure 2: T-shape structure 

are used to solve the substructured problems, respectively used together with 
Neumann-Neumann and Dirichlet preconditioners. Beside, the convergence 
criterion of the solver, which stands here for the interface traction gap (resp. 
displacement gap) in the primal (resp. dual) approach, is set to 10~ 6 . 

On each case, in addition to the new parallel error estimator e^ 11 , we 
compute the standard sequential e^ and the true error obtained using a 
reference field u ex computed on a very fine mesh: 



e CR = e CR(n)(Uh,(?h) 




6.1. Quality of the parallel error estimator 

We first study the quality of the parallel error estimator for com- 
putations when convergence of the domain decomposition solver is reached. 
As said earlier, the proposed technique does not lead to the same statically 
admissible field because of the special treatment of the interface traction 
(COI]). Our estimator might then be sensitive to the substructuring, we thus 
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compare the estimations obtained with meshes of characteristic size h and 
decomposition into N sd subdomains. Results are given in Figure Eland Table 

m 



Exact error and estimate as function of mesh size 
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Figure 3: Convergence of Error e/j and Estimators e^S and e^ 11 (for various N s d) vs. 
element size h 

We observe that: 

• The results obtained by FETI and BDD can not be distinguished 
(which is why only FETI results are given in Tabled]). 

• e cR 1 barely depends on the substructuring; the results are quite similar 
whether they are conducted on a single domain ( "sequential" curve) or 
on N s d subdomains. Only a slight rise of the estimation can be observed 
when the number of interface degrees of freedom is not small compared 
to the number of internal degrees of freedom, which is logical since 
the description of interface traction fields is coarser in parallel than in 
sequential. 

As a conclusion, the parallel error estimator e^ 11 enables to recover the 
same efficiency factor as the standard sequential one, while the CPU-time is 
divided by N sd . 
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h 


L/2 


L/4 

/ 


L/8 


L/16 


L/32 


# dofs 


146 


514 


1922 


7426 


29186 


e h 


0.2347 


0.1493 


0.0937 


0.0597 


0.0386 


e CR 


0.5712 


0.4035 


0.2662 


0.1769 


0.1151 


N sd 






„ddm 
e CR 






2 


0.5657 


0.4021 


0.2648 


0.1747 


0.1151 


4 


0.5768 


0.4007 


0.2648 


0.1747 


0.1151 


8 


0.5546 


0.4007 


0.2676 


0.1747 


0.1165 


16 






0.2690 


0.1761 


0.1165 


32 






0.2787 


0.1789 


0.1178 



Tabic 1: Error et and Estimators and e^ 11 (for various N s d) vs. element size h 



6.2. Convergence of the parallel estimator along DD-solver iterations 

Previous results enabled to analyse the quality of the parallel estimator 
when interface quantities had converged. A new feature associated to the use 
of an iterative solver for the domain decomposition (DD) problem is that the 
discretization error estimation can be conducted before DD convergence is 
reached, that is in presence of displacement or traction discontinuity at the 
interface as explained in Sections H] and |5] 

We then compute the parallel error estimator e^f at each iteration of the 
DD solver. Convergence curves of e^" during the FETI and BDD iterations 
are shown on Figure HJ Parallel error estimator is plotted as a function of 
the FETI (resp. BDD) residual, defined (for Iteration n) as the normalized 
displacement (resp. traction) gap at the interface: 

" Mr " Mr (30) 

Classical stopping criterion for the of convergence of DD solver is this residual 
being below 10~ 6 . Because of the similarity between the curves, the only 
shown cases correspond to h = L/8 and h = L/16. 

The curves show a rapid convergence of the parallel error estimator along 
iterations of the solver, so that can be considered as converged when 
FETI residual reaches an order of magnitude of 5.10~ 3 or BDD residual 
reaches 5.10" 1 , which corresponds to at most 5 iterations whereas the solver 
convergence is achieved in 10 to 20 iterations. 
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Figure 4: Convergence of estimator vs. DD residual 



Actually, e^ 11 is driven by both the discretization error and the conver- 
gence of the solver (interface error). The "L"-shaped curves show that the 
impact of residual of the DD solver is preponderant only at the first iterations 
(when interface fields are very poorly estimated), after e^ 11 stagnates at a 
value very close to e^ which is only associated to the discretization error. 

Then, it seems possible to stop the iterations of the solver far before 
convergence while still obtaining an accurate global estimate for the dis- 
cretization error. 

Figures [5] and [6] show maps of the elementary contributions e^ 11 ^ to the 
parallel error estimator e^ 11 at different steps of the convergence, for N S d = 8 
with h e = L/2 or h e = L/8. At the first iterations it can be seen that the 
estimator highlights both discretization errors (around the re-entrant angle) 
and lack of convergence of the solver (along the interfaces), whereas very 
quickly the solver (that is the interfaces) does not contribute any more to 
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the estimator. 

The various examples show that the convergence of the global estimator is 
due to the convergence of elementary Contributions e^[ n ' E , which means that 
when willing to carry out remeshing procedures, the maps obtained after few 
iterations of the solver are sufficient to define correct refinement instructions. 

7. Conclusions 

In this paper, we presented a new approach to handle robust model veri- 
fication based on constitutive relation error in a domain decomposition con- 
text. 

The method relies on the construction of fields that are kinematically 
and statically admissible on the whole structure. We showed that a fully 
parallel construction is possible even when starting from fields which do not 
satisfy interface conditions. The construction is a three-step procedure: first 
displacement and traction nodal fields are built-up so that discrete admis- 
sibility conditions are satisfied, second continuous admissible traction fields 
are deduced, third these fields are used as input by any classical recovery 
procedure. The first step is implicitly done when good preconditioners are 
employed within the domain decomposition methods and the second step 
corresponds to the inversion of small and sparse "mass" matrices. 

Our first results show that not only the estimation error does not suffer 
from the approximation that are made at the interface in order to achieve full 
parallelism, but that even roughly estimated interface fields enable to obtain 
a good estimation of the discretization error and correct maps of elementary 
contributions which are required by mesh adaptation procedures. Thus not 
only the computational cost are divided by the number of processors but 
the prior obtainment of the finite element solution can be accelerated since a 
coarse solution is sufficient (which corresponded to 3 to 5 times less iterations 
in our case). 

Future studies will deal with parallel mesh adaptation. 
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Figure 5: Maps of e^ n ' E for h = L/2 and N s d = 8 at various iterations 
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Figure 6: Maps of e^p 111 ' 51 for h = L/8 and N s d = 8 at various iterations 



24 



Id) 2W 3 (D 




m — 3(2) 

(a) Subdomains 





(b) Local interface 




(c) Primal interface 



(d) Dual interface 



A (i) 
A(D 



?0 








t (2) 



A (2) 



A (2) 




t (3) 



A (3) 



A (3) 



'10 
1 
\0 



/0 

1 


/ 



-1 




V o 



1 





1 Oy 












oj 



Figure 7: Local numberings, interface numberings, trace and assembly operators 



25 



Fast estimation of discretization error for FE problems 
solved by domain decomposition 



bo 



A.Parret-Freaud a , C. Rey*' a , P. Gosselet a , F. Feyel b 

a LMT-Cachan, ENS Cachan/CNRS/UPMC/PRES UniverSud, 61 av. du president 

Wilson, 94235 Cachan cedex, France 
b ONER A, DMSM/CEMN, 29 avenue de la division Leclerc, BP72, F92322 Chatillon 

cedex, France 



Abstract 

This paper presents a strategy for a posteriori error estimation for substruc- 
tured problems solved by non-overlapping domain decomposition methods. 
We focus on global estimates of the discretization error obtained through the 
error in constitutive relation for linear mechanical problems. Our method al- 
lows to compute error estimate in a fully parallel way for both primal (BDD) 
and dual (FETI) approaches of non-overlapping domain decomposition what- 
ever the state (converged or not) of the associated iterative solver. Results 
obtained on an academic problem show that the strategy we propose is effi- 
cient in the sense that correct estimation is obtained with fully parallel com- 
putations; they also indicate that the estimation of the discretization error 
reaches sufficient precision in very few iterations of the domain decomposi- 
tion solver, which enables to consider highly effective adaptive computational 
strategies. 

Key words: verification; error in constitutive relation; non overlapping 
domain decomposition; FETI; BDD. 



1. Introduction 

The setting-up of robust numerical methods to solve complex systems of 
partial differential equations has become a key issue in applied mathematics 
and engineering, driven by the increasing use of numerical simulation in both 
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research and industry. Among the latter, virtual testing has become a short 
term aim, with the objective to replace expensive experimental studies and 
validations by numerical simulations, even in order to certify large structures 
as planes and bridges. 

Thus, one key point of the numerical methods to develop is the verification 
of computations which enables to warranty that the computed solution is 
sufficiently close to the original continuum mechanics model. This topic of 
numerical analysis has been the subject of many studies for the last decades. 
Three main classes of error estimator have been developed, based either on 
equilibrium residuals [? ], flux projection [? ] or error in constitutive law [? 
]. An overview of those various methods can be found in [? ]. 

Another key point of numerical methods is their ability to quickly provide 
solutions to large (nonlinear) systems. The most classical answer to this issue 
is to use domain decomposition methods in order to take advantage of the 
parallel hardware architecture of recent clusters and grids. In engineering, 
non-overlapping domain decomposition methods are mostly employed, such 
as the well known FETI [? j or BDD [? ]. An overview of the main 
approaches related to non-overlapping domain decomposition can be found 
in [? ]. 

We aim to provide fully integrated adaptive strategies to compute large 
structural mechanics problems with certified quality. To do that, our cur- 
rent approach is to explore some ways of making bidirectional interactions 
between domain decomposition and a posteriori error estimation. Our de- 
velopments are based both on the error in constitutive relation to measure 
the quality of our results and to forecast mesh refinement, and on a generic 
vision of non-overlapping domain decomposition methods which enables to 
do high-performance computing. 

This paper focuses on the estimation of the global error in constitutive re- 
lation in order (among others) to study how it is influenced by the error in the 
convergence of the domain decomposition solver which is linked to the non- 
satisfaction of interface equations (continuity of displacements and balance 
of forces). To do so we propose a strategy to build, in parallel and during the 
iterations, displacement and stress fields which are kinematically admissible 
(KA) and statically admissible (SA) on the whole structure. We face two 
main difficulties. First, since before convergence interface fields do not pos- 
sess the classical properties of discretized fields (continuity of displacements 
and weak equilibrium), the recovery of admissible displacements and stresses 
requires some preprocessing. Second, the computation of statically admissi- 
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ble fields being an operation which can not be conducted independently on 
each element (in some methods it can even be a large bandwidth operation), 
classical recovery methods [??????] would require inter-sub domain 
communications. 

Our generic method to build continuous displacement and balanced trac- 
tion fields for both primal and dual approaches of non-overlapping domain 
decomposition is presented through this paper. It will be shown that the 
properties of the preconditioners involved in domain decomposition solvers 
make this reconstruction costless, and that an error estimator can then be 
computed in a fully parallel way. 

This paper is organized as follows. Section [5]recalls the general framework 
related to our upcoming developments, mainly the estimation of the error in 
constitutive equation and the use of domain decomposition method. Section 
[3] shows how the problem of error estimation in a substructured context can 
be brought back to the computation of nodal displacement and traction fields 
which are admissible in a discrete sense. Sections H] and |5] describes how to 
obtain these fields without inter-sub domains exchanges when using classical 
primal (BDD) and dual (FETI) domain decomposition methods with good 
preconditioners. Section [6] presents numerical assessments, first to validate 
the parallel recovery procedure, then to prove that a good estimation can 
be obtained far earlier than the solver converged (in the sense of domain 
decomposition iterative solver). Finally, Section [7] concludes this paper. 

2. Framework of the study 

2.1. Reference mechanical problem 

Let us consider the static equilibrium of a structure which occupies the 
open domain Q C W 1 and which is submitted to given body forces /, to given 
traction forces g on dfQ and to given displacements uq on the complementary 
part d u fl 7^ 0. We assume the structure undergoes small perturbations and 
that the material is linear elastic, characterized by the Hooke's tensor EL 
Let u be the unknown displacement field, e(u) the symmetric part of the 
gradient, a the Cauchy stress tensor. 

Let u) C Q be an open subset of Q, dfU = dui fl dffl, d u u = du n d u Q and 
r = du \ (d u u U dfu) (see Figured]). We introduce two affine subspaces and 
one positive form: 
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Figure 1: Domain Q,, subdomain lo and boundaries 

• Subspace of kinematically admissible fields 

KA(w) = jw G (H 1 (a;)) <i , tr(-u) = u on <9 M wj 

where tr is the trace operator. 

• Subspace of statically admissible fields 

SA(w) = |r G (L 2 (cu)) dxd ,r symmetric, Vw* G KA 00 (w), 

/ r : e(u*)du = / f.u*du + / g.u*dS 

where KA°V) = |u G (H 1 ^))^ , tr(u) = on d u co U r| 

• Measure of the non-verification of the constitutive equation [? ] 

ecR(«)(u, <r) = ||<t - H : e (u) || H -i, w 



where llxlL-i 



^ (x : H^ 1 : x) doo 



The mechanical problem set on Vt can be formulated as: 

Find {u ex ,a ex ) G KA(fi) x SA(Q) such that e C R(n)(u ex , <J ex ) = 

4 



2.2. Finite element approximation for the global problem 

Let Qh be a tessellation of Cl to which we associate a finite dimensional 
subspace KA^(f2) of KA(fi). The classical finite element displacement ap- 
proximation consists in searching 

u h G KA h (n) 
a h = H : e{u h ) 

a h : e(u* h )dQ = [ f.u* h dQ + [ g.u* h dS, V< G KA° h °(Q) 



n Jn Jd f n 

After introducing the d x A^o/ matrix ip h of shape functions which form 
a basis of KA/j(f2) and the vector of nodal unknowns u (of size A^o/, number 
of degrees of freedom) so that Uh = ¥>h u > the classical finite element method 
leads to the well-known linear system: 

Ku = f (5) 

where K is the (symmetric positive definite) stiffness matrix of domain fl^ 
and f is the vector of generalized forces. 

2.3. A posteriori error estimator 

The finite element approximation (uh, ah) satisfies Uh G KA(Q) and 
e CR(a)(uh, o~h) = but ah ^ SA(fi). The error in constitutive relation 
consists in deducing from (uh, ah) an admissible displacement-stress pair 
(uh,ah) G KA(Q) x SA(fi) in order to measure the residual on the consti- 
tutive equation ([3]) ecR(o)(«/i, a h ) ^ 0. Using the well-known Prager-Synge 
theorem it can be proved that 

\\e(u ex ) -e(wfc)llH,n+ Ike* - ^lln-i.n = { e CR(n){u h ,Vh)) 2 

Hence, the evaluation of the error in constitutive relation ecR(n)(uh,&h) for 
any admissible pair (uh, a~h) provides a guaranteed upper bound of the global 
error 

\\e{u ex ) -e(u A )||H,n < e CR{n) (uh,a h ) (6) 

KA/ l (f2) being a subspace of KA(f2), the construction of an admissible dis- 
placement field Uh is straightforward since it can be taken equal to Uh- On the 
other hand, as ah is not statically admissible, the construction of an admis- 
sible stress field dh G SA(f2) is a crucial point which has already been widely 
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studied in the literature. A first solution is to use a dual formulation of the 
reference problem [? ] to compute a from scratch. Unfortunately building a 
subspace of SA(fi) is a complex task and most people prefer to post-process 
a statically admissible field from Field ah obtained by a displacement for- 
mulation. Classical methods are the element equilibration techniques [? ? 
], which have been improved by the use of the concept of partition of unity 
which lead to [? ? ? ] and the flux- free method [? ]. In most cases they 
involve the computation of efforts on "star-patches" which are the set of el- 
ements sharing one node, for each node of the mesh. Though rather simple 
these computations are in great number and thus expensive. 

In the following, we note by T h the algorithm which has been chosen to 
build an admissible stress field a^. Whatever the choice, the algorithm takes 
as input not only the finite element stress field ah but also the continuous 
representation of the imposed forces (f,g). 

ok = ?h{?h,f,g) e SA(fi) 

The algorithm we have used for our applications is the one proposed in [? ] 
using a three degrees higher polynomial basis when solving the local problems 
on elements [? ]. 

2.4- Substructured formulation 

Let us consider a decomposition of domain Vt in open subsets (^^)i< s <jv sd 
(N sd is the number of subdomains) so that £1^ fl = for s ^ s' and 
0, = U s f^ s ). Let u u = (u^)si we define the global assembling operator A: 

u = A(u a ) & u ln{s) = u {s) (7) 

In order to reformulate the mechanical problem on the substructured con- 
figuration, we need to specify the conditions that should be satisfied at the 
boundary between subdomains r^ s '- ) = dQ^ fl dil^ s '\ We have the funda- 
mental properties: 

E KA ^ * { tr(«W) = tr(«<W) on V™, V( S , s>) (8) 
A{<P) e SA(fl) * | ffW _ nW + \ {s { n(sl) = Q on r w V(S) g/) (9) 
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In other words, in order to be admissible on the whole domain fi, not only 
fields need to be admissible in a local sense (independently on each f2^), but 
they also need to satisfy interface conditions, namely displacements continu- 
ity and tractions balance (action-reaction principle). 

2.5. Finite element approximation for the sub structured problem 

We assume that the tessellation of Cl and the substructuring are conform- 
ing so that (i) each element only belongs to one subdomain and (ii) nodes are 
matching on the interfaces. Each degree of freedom is either located inside 
a subdomain (subscript i) or on its boundary = U S >T( SS > (subscript b) 
where it is shared with at least one neighboring subdomain. Let A^ be the 
vector of unknown efforts imposed on the interface of subdomain fl^ by its 
neighbors. The finite element problem (jSJ) can be written highlighting the 
contributions of subdomains: 



Ea^a^ = o 



Ws, K^uW = f« + t« J A« with { £ a(s)u ( s) = q (10) 

s 

where is the discrete trace operator (u£ = t^u^) and where 
and are assembling operators so that A^ s ^ enables to formulate the 
mechanical equilibrium of interfaces (Q and A^ enables to formulate the 
continuity of displacements dSJ (in the case of two subdomains, we have 
Es Ais) *b = A[ 1} + A[ 2) = and £sA (s) i4 s) = u < 1} - u< 2) = 0, see Fig. [7J 
for less trivial example and [? ] for more an extensive description of all 
operators). One fundamental property of assembling operators is their or- 
thogonality: 

^aWaw t = o (11) 

s 

Note that the equilibrium of subdomain also writes: 

or in an equivalent condensed form: 

S«iiW=bW+AM (13) 
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with 

where SW is the Schur complement and b( s ) is the condensed right-hand side. 



3. A posteriori error estimator in substructured context 

The key point for the efficient evaluation of the error in constitutive rela- 
tion in a substructured context (without overlapping) is to define admissible 
pairs (u^\a^) € KA(f^) x SA(f^) on each subdomain so that the associ- 
ated assembled pair is admissible for the reference problem (A(v^), A(a^)) G 
KA(Q) x SA(fi). Due to the absence of overlap, the additive structure of the 
associated error in constitutive relation leads to a fully parallel evaluation of 
the a posteriori error estimator: 

(e C R(n){A(^),A(dh))) 2 = Yl ( e CR(f^)) ^ s) )) 

s 

The application of a classical recovery strategy to compute admissible 
fields raises two difficulties in a substructured context. First, the star-patches 
can not be employed on the boundary nodes without assuming communica- 
tion between subdomains. Though these exchanges would remain limited, 
we propose an alternate strategy to achieve full parallelism without impair- 
ing the properties of the error in constitutive relation. Second, in order to 
solve the substructured problem (fit)]) parallel strategies consist in using iter- 
ative solvers which are based on the loosening of at least one of the interface 
conditions which is only verified (up to a certain precision) once the solver 
converged. Thus recovering strategies need to be adapted so that the local 
fields (uj*\dfr) satisfy the interface conditions. 

The aim of this section is to prove that the determination of the admissi- 
ble pair (A(u^), A(&h)) can be brought back to the determination of nodal 

interface fields (u£ s \A^) s which satisfy specific interface conditions. The 
construction of these nodal fields depends on the chosen domain decomposi- 
tion strategy and is discussed in the following sections. 
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3.1. Kinematically admissible fields 

In order to ensure interface Condition (JS} when building G KA(f2^) 
so that A(v^) G KA(f2), we introduce continuous interface displacement 

A (S) 

fields from which we shall deduce internal displacement fields: 



= fig? , v(«y) 

?y (s) - ?y (s) Vs 
fc|r( ss ') ~~ 6^ ' 

Since discretizations are matching on the interface, the first condition can 
directly be imposed on finite element nodal quantities: 

u< s) = uj s '\ V( S , S ') 

In order to deduce the internal fields, one finite element problem is solved 
independently on each subdomain with imposed Dirichlet conditions on the 
interface: 



a!" = kF 



U 6 



tt = -4(^) G KA(fi) 

S 1 .^. Statically admissible fields 

In order to ensure interface Condition ([5]) when building a^' G SA(f2^) 
so that A(Oh) G SA(O), we introduce for each subdomain the continuous 
balanced interface traction fields defined on which satisfy: 

^ ] -n {s) = on rW 

*£> + = o on r(-') (15) 

/ f.pdtt+ [ g('\pdS+ [ F$.pdS = Vp G RKA°(fiW) 

where RKA°(f2( s )) is the set of rigid body motions which are compatible with 
Dirichlet conditions imposed on d u Cl^: 

RKA°(fi (s) ) = {pe H^fiW), p = on d u ^ s) , e(p) = 0, } 
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The last condition of (I15p is the translation of Fredholm's alternative in 
order to ensure the well-posedness of the static problem on domain fiw. To 
build these traction fields in a simple way, we associate them with the finite 



element nodal reaction field X b : 

/ ,^-^\ r ^dS = \% (16) 
where j denote a node of the interface, ipj its associated shape function and 

~(s) ~(s) 

\ b j the corresponding nodal component of X b . This equation then imposes 

-~(s) ^(s) 

that the discrete field X b and the continuous field F bh develop the same 
virtual work in any finite element displacement field. The conditions on F bh 

have these discrete counterparts on X b : 

£ A( ^ S) = ° 

(17) 

R (s) T (t(°) T \ ( b s) + f M) = o 

where is a basis of ker(K^^). As said earlier, the first equation corre- 
sponds to the equilibrium between subdomains. The second equation corre- 
sponds to the balance of the subdomain with respect to virtual rigid body 
motions (since this kind of displacement field is exactly represented in the 
finite element approximation, the discrete condition is equivalent to the con- 
tinuous one). 

As a first approach, we define F b ^ as: 

F« = ,F« (18) 

"^(b) ~^(s) (s) 

where F b is the vector of nodal values of F bh and <Ph) T [ a ) refers to the vector 

of the trace on of finite element shape functions. Vector F b is then 
obtained by the inversion of the (small) "mass" matrix of the interface of 
each subdomain. In the following, we denote by Qh the previous procedure 
which associates a continuous balanced interface force F b ^ to a balanced 

~(s) 

nodal interfaces forces X b : 

F b (s h ] = e h & s) ) (is) 
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The traction field F b ^ allows to satisfy the interface conditions associated 
to the static admissibility. The next step is to build internal finite element 

~(s) 

stress fields which match the associated nodal boundary field X b . This is 
done by solving one finite element problem on each subdomain with imposed 
Neumann conditions on the interface. 

u« = K« + (f W + t^A^) 

^=J-,(H: £ (^uW),/«,{^),^(Ah}) ^ 
a h = Affi) e SA(fl) 

The use of the pseudo-inverse is due to the potential lack of Dirichlet 

boundary conditions on the substructure. Displacement field is defined 
up to a rigid body motion which needs not to be determined since only the 
symmetric gradient of the associated displacement field is required. 

It has to be noted that the fully parallel procedure Gh proposed above 
leads to a different admissible traction field as would have been obtained 
using standard patch-technique [? ] (referred in the sequel as the sequential 
approach). Thus the use of Gh implies that the parallel error estimation 
is different from the standard sequential one even when discrete interface 
conditions are satisfied. For now there are no theoretical results on the 
quality of the resulting fields, examples (as given in Section [6]) show that 
sequential estimator and parallel estimator (when interface conditions have 
sufficiently converged, which happens very quickly) can not be distinguished. 

4. Recovery of admissible fields in BDD 

In the Balancing domain decomposition [? ? ], a unique interface dis- 
placement unknown U(, is introduced so that continuity is always insured: 

u[ s) = A^ T u b =► A (s) ui s) = (21) 

s 

Other quantities can be deduced from u b and equations (I12II13P : 
\ { b s) = sWu6-b« 
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The BDD solver consists in iteratively finding the interface displacement U;, 
which insure global equilibrium (Y^ s A^A^ = 0), 



o = A(s)x l s) = (x> (s)s(s)A(s)T ) u * - ( E A(s)b(s) 

S \ S / \ S / 



(23) 



4-1. Recovery of KA fields 

In the BDD solver, kinematic interface conditions are satisfied anytime 
and using = u;, enables to build uf 3 so that Uh = A (v^) e KA(fi). Note 
that all associated computations are realized during the standard resolution 
process so that no extra operation is required. 

4-2. Recovery of SA fields 

For a given interface displacement we note: 

M = £ a«a<'> = A(s) ( s(s) ^ - b(s) ) 

s s 

Obviously |JAfoJ| is zero if and only if u b is the solution to (123"]). We then 
define: 

A^A^- A« T M (24) 

where (A^) s are scaled assembling operators so that ^ s A^) A^) = I. 
The multiplicity scaling is a typical example of such operator AW; 

T 1 

which, in the case of two subdomains, gives A^ [fA&]] = — ff Af,7l • I* 1 t ne 
case of heterogeneous structures, other scaled assembly operators which take 
the heterogeneity into account are used [? ? ? ]. 

~(s) 

It is clear that by definition, X b is a balanced nodal reaction field: 

s 

~(s) 

In order to prove that A b also satisfies Fredholm's alternative, we note that 
since is a basis of ker(K^) and is invertible, we have S^R^ = 
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and R 4 = — K)* J K^'R^ . The condition then writes in an equivalent 
condensed form: 

^ ,T (Ar +b ">)=o 

R< S)T (s'-'uj - b<»» + A«'> T M + b<*>) = 

Using the symmetry of (inherited from the symmetry of K^) to nullify 
H b S^ sj , the condition writes: 

(AWR( S) ) T JA 6 H = (25) 

which is exactly the balancing condition [? ] of the iterative BDD solver: the 
residual of the BDD iterative solver f[A 6 "fl = fe s A^aJ s) ) ([23]) has to be 
orthogonal to all local weighted rigid body motions so that preconditioning 
step is well posed. 

Then we have constructed a pair of interface nodal Vectors (u&, X b ) which 
satisfy all required conditions to build admissible fields. 

Note that all the involved operations are already realized during clas- 
sical steps of the primal domain decomposition approach with a Neumann- 
Neumann preconditioner and the associated coarse problem, so that all finite 
element quantities (even the internal ones) are available at no cost; the only 
extra operations are due to the use of Algorithms Qh (to compute F bh ) and 
Th (to compute ah). 

5. Recovery of admissible fields in FETI 

In the Finite Element Tearing and Interconnecting domain decomposition 
[? ], a unique interface effort unknown At, is introduced so that interface 
equilibrium is always insured: 

\ { b s) = A« T A b A(s)A i S) = ( 26 ) 

s 

Displacements can be deduced from X b if it satisfies Fredholm's alternative 
on each substructure: 

U W = (f« + t^ T A^ T \ b ) + R«aM 

V (27) 
= RW r (f « + t^ T A^ T \ b ) 
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where ct^ is the unknown magnitude of rigid body motions. The FETI 
solver consists in iteratively finding an interface effort A&, under the previous 
constraint, which insures the continuity of interface displacement: 

o = J2 A (s)u i s) = fe a^k^\^ t a^A \ b 

+ ^A^K^A + ^AW'WwA (28) 

5.1. Recovery of SA fields 

In the FETI solver, the nodal interface fields A^ = A^ T A& are by 
construction always balanced at the interface (1261) and associated to well- 
posed discrete Neumann problems on each substructure (l2Tj) . Hence, we 



can directly set a[ ^ = A^ and apply algorithms Gh an d J~h to compute 
d ( h s) e SA(Q^) with a h = Affi) e SA(fl). 

5.2. Recovery of KA fields 

For a given balanced nodal interface traction A&, we introduce, in agree- 
ment with ( 1271) . the gap of the interface displacement : 

s 

and we define 

~( s ) 00 aW T m II 

u b = u l - A M 

~ (s) / \ ~ (s) T 

where (A ) a are scaled assembling operators so that Y]„ A^ s) A = I. 

~ (s) 

Similarly to the BDD typical example of such operator A is the 

multiplicity scaling: 



-i 



A (s) = A^ T [J2A U) A^ T ) (29) 



~ (s) T 1 

Note that in the case of two subdomains, we have: A [[ubJJ = — [|u&J|- 

The connection between FETI and BDD scaling operators (even in the het- 
erogeneous case) is given in [? ]. 
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It is clear that by construction 

s 

Hence nodal interface displacement can be used to deduce an admissible 
displacement field uj^ so that Uh = A (u^) G KA(fi). 

Then we have constructed a pair of interface nodal Vectors (uf,, A&) which 
satisfies all required conditions to build admissible fields. Note that all the 
involved operations are already realized during classical steps of the dual do- 
main decomposition approach (with built-in coarse problem) with Diri chief's 
preconditioner, so that all finite element quantities (even the internal ones) 
are available at no cost: the quantity [[u^JJ is directly available during the 
classical solution procedure (without computing any afity which is based on 
an initialization/projection algorithm [? ], and the displacement field 
can be defined up to an element of the kernel (a rigid body motion) since 
only its symmetric gradient is used during the computation of the error. The 
only extra operations are due to the use of algorithms Gh (to compute F^h) 
and Th (to compute ah). 



6. Numerical assessment 

In order to assess the performance of our parallel error estimator, we 



consider the 2D toy problem of the T-shape structure of Figure 2(a) which has 
been used in other papers like [? ] . Plane stresses are assumed. The material 
behavior is isotropic, linear and elastic, with Young modulus E = 2000 MPa 
and Poisson's ratio v = 0.3. The structure is clamped on its basis (whose 
length is denoted L) and it is submitted to traction and shear on its upper- 
right side, while all the remaining boundaries are traction-free. 

Several regular meshes have been generated, constituted by triangular 
elements of characteristic size h = - with m = 2,4,8,16,32. For each 

m 

mesh, a sequential (mono-domain) computation is driven, followed by domain 
decomposition computations obtained by an automatic splitting of the mesh 
in an increasing number N s d of subdomains ( N s d = 2, 4, 8 when m ^ 4 and 



N s d = 2,4,8, 16,32 when m ^ 8). Figure 2(b) shows such a decomposition 
for N s d = 8 and m — 8. 

All the computations are driven in the ZeBuLoN finite element code [? ], 
using elements of polynomial degree p = 1. Both BDD and FETI algorithms 
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(a) Finite element problem (b) Substructuring 

(h= §) (h= ±,N sd = 8) 



Figure 2: T-shape structure 

are used to solve the substructured problems, respectively used together with 
Neumann-Neumann and Dirichlet preconditioners. Beside, the convergence 
criterion of the solver, which stands here for the interface traction gap (resp. 
displacement gap) in the primal (resp. dual) approach, is set to 10~ 6 . 

On each case, in addition to the new parallel error estimator e^ 11 , we 
compute the standard sequential e^ and the true error obtained using a 
reference field u ex computed on a very fine mesh: 



e CR = e CR(n)(Uh,(?h) 




6.1. Quality of the parallel error estimator 

We first study the quality of the parallel error estimator for com- 
putations when convergence of the domain decomposition solver is reached. 
As said earlier, the proposed technique does not lead to the same statically 
admissible field because of the special treatment of the interface traction 
(COI]). Our estimator might then be sensitive to the substructuring, we thus 
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compare the estimations obtained with meshes of characteristic size h and 
decomposition into N sd subdomains. Results are given in Figure Eland Table 

m 



Exact error and estimate as function of mesh size 
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Figure 3: Convergence of Error e/j and Estimators e^S and e^ 11 (for various N s d) vs. 
element size h 

We observe that: 

• The results obtained by FETI and BDD can not be distinguished 
(which is why only FETI results are given in Tabled]). 

• e cR 1 barely depends on the substructuring; the results are quite similar 
whether they are conducted on a single domain ( "sequential" curve) or 
on N s d subdomains. Only a slight rise of the estimation can be observed 
when the number of interface degrees of freedom is not small compared 
to the number of internal degrees of freedom, which is logical since 
the description of interface traction fields is coarser in parallel than in 
sequential. 

As a conclusion, the parallel error estimator e^ 11 enables to recover the 
same efficiency factor as the standard sequential one, while the CPU-time is 
divided by N sd . 



17 



h 


L/2 


L/4 

/ 


L/8 


L/16 


L/32 


# dofs 


146 


514 


1922 


7426 


29186 


e h 


0.2347 


0.1493 


0.0937 


0.0597 


0.0386 


e CR 


0.5712 


0.4035 


0.2662 


0.1769 


0.1151 


N sd 






„ddm 
e CR 






2 


0.5657 


0.4021 


0.2648 


0.1747 


0.1151 


4 


0.5768 


0.4007 


0.2648 


0.1747 


0.1151 


8 


0.5546 


0.4007 


0.2676 


0.1747 


0.1165 


16 






0.2690 


0.1761 


0.1165 


32 






0.2787 


0.1789 


0.1178 



Tabic 1: Error et and Estimators and e^ 11 (for various N s d) vs. element size h 



6.2. Convergence of the parallel estimator along DD-solver iterations 

Previous results enabled to analyse the quality of the parallel estimator 
when interface quantities had converged. A new feature associated to the use 
of an iterative solver for the domain decomposition (DD) problem is that the 
discretization error estimation can be conducted before DD convergence is 
reached, that is in presence of displacement or traction discontinuity at the 
interface as explained in Sections H] and |5] 

We then compute the error estimator at each iteration of the DD solver. 
Convergence curves of e^" during the FETI and BDD iterations are shown 
on Figure HJ Error estimator is plotted as a function of the FETI (resp. 
BDD) residual, defined (for Iteration n) as the normalized displacement (resp. 
traction) gap at the interface: 

" Mr " Mr (30) 

Classical stopping criterion for the of convergence of DD solver is this residual 
being below 10~ 6 . Because of the similarity between the curves, the only 
shown cases correspond to h = L/8 and h = L/16. 

The curves show a rapid convergence of the parallel error estimator along 
iterations of the solver, so that can be considered as converged when 
FETI residual reaches an order of magnitude of 5.10~ 3 or BDD residual 
reaches 5.10" 1 , which corresponds to at most 5 iterations whereas the solver 
convergence is achieved in 10 to 20 iterations. 
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Figure 4: Convergence of estimator vs. DD residual 



Actually, e^ 11 is driven by both the discretization error and the conver- 
gence of the solver. The "L" -shaped curves show that the impact of residual 
of the DD solver is preponderant only at the first iterations (when interface 
fields are very poorly estimated), after e^ 11 stagnates at a value very close 
to e^pl which is only associated to the discretization error. 

Then, it seems possible to stop the iterations of the solver far before 
convergence while still obtaining an accurate global estimate for the dis- 
cretization error. 

Figures [5] and [6] show maps of the elementary contributions e^ 11 ^ to the 
parallel error estimator e^ 11 at different steps of the convergence, for N s d = 8 
with h e = L/2 or h e = L/8. At the first iterations it can be seen that the 
estimator highlight both discretization errors (around the re-entrant angle) 
and lack of convergence of the solver (along the interfaces), whereas very 
quickly the solver (that is the interfaces) does not contribute any more to 
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the estimator. 

The various examples show that the convergence of the global estimator is 
due to the convergence of elementary Contributions e^[ n ' E , which means that 
when willing to carry out remeshing procedures, the maps obtained after few 
iterations of the solver are sufficient to define correct refinement instructions. 

7. Conclusions 

In this paper, we presented a new approach to handle robust model veri- 
fication based on constitutive relation error in a domain decomposition con- 
text. 

The method relies on the construction of fields that are kinematically 
and statically admissible on the whole structure. We showed that a fully 
parallel construction is possible even when starting from fields which do not 
satisfy interface conditions. The construction is a three-step procedure: first 
displacement and traction nodal fields are built-up so that discrete admis- 
sibility conditions are satisfied, second continuous admissible traction fields 
are deduced, third these fields are used as input by any classical recovery 
procedure. The first step is implicitly done when good preconditioners are 
employed within the domain decomposition methods and the second step 
corresponds to the inversion of small and sparse "mass" matrices. 

Our first results show that not only the estimation error does not suffer 
from the approximation that are made at the interface in order to achieve full 
parallelism, but that even roughly estimated interface fields enable to obtain 
a good estimation of the discretization error and correct maps of elementary 
contributions which are required by mesh adaptation procedures. Thus not 
only the computational cost are divided by the number of processors but 
the prior obtainment of the finite element solution can be accelerated since a 
coarse solution is sufficient (which corresponded to 3 to 5 times less iterations 
in our case). 

Future studies will deal with parallel mesh adaptation. 
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Decomposition Reference map 
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Iteration 4 
BDD solver 



Iteration 5 




Iteration 1 



Iteration 4 
FETI solver 
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Figure 5: Maps of e^ n ' E for h = L/2 and N s d = 8 at various iterations 
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Decomposition Reference map 
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Iteration 1 



Iteration 4 
BDD solver 
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Iteration 1 



Iteration 4 
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Iteration 5 



Figure 6: Maps of e^p 111 ' 51 for h = L/8 and N s d = 8 at various iterations 
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Figure 7: Local numberings, interface numberings, trace and assembly operators 
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